Time-resolved extinction rates of stochastic populations 
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Extinction of a long-lived isolated stochastic population can be described as an exponentially slow 
decay of quasi-stationary probability distribution of the population size. We address extinction of 
a population in a two-population system in the case when the population turnover - renewal and 
removal - is much slower than all other processes. In this case there is a time scale separation in the 
system which enables one to introduce a short-time quasi-stationary extinction rate Wi and a long- 
time quasi-stationary extinction rate Wi, and develop a time-dependent theory of the transition 
between the two rates. It is shown that Wi and Wi coincide with the extinction rates when the 
population turnover is absent, and present but very slow, respectively. The exponentially large 
disparity between the two rates reflects fragility of the extinction rate in the population dynamics 
without turnover. 

PACS numbers: 02.50.Ga, 87.23. Cc 
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I. INTRODUCTION 

An isolated stochastic population (of molecules, bacte- 
ria, animals, parasites inhabiting a community of hosts, 
etc.) ultimately goes extinct. The ultimate extinction is 
driven, even in the absence of unfavorable environmen- 
tal variations, by large demographic fluctuations: chains 
of random events when population losses dominate over 
gains. The risk of population extinction is a major is- 
sue in viability of small populations in the nature 0, Q , 
whereas extinction of an endemic disease from a commu- 
nity |3[ is a desirable development. 

Calculating the extinction rate of a stochastic popu- 
lation is a challenging problem. Here one needs to eval- 
uate the (very low) probability of a large fluctuation in 
a stochastic system which is far from equilibrium and 
therefore defies many standard assumptions and methods 
of statistical mechanics. In spite of this difficulty a sig- 
nificant progress has been achieved in this type of prob- 
lems via the use of a WKB (Wentzel-Kramers-Brillouin) 
approximation borrowed from quantum mechanics (or, 
more generally, wave mechanics) and adapted to the dis- 
sipative, non-Hermitian stochastic Markov processes Q. 
The WKB approximation employs, as a large parame- 
ter, the typical population size in the metastable quasi- 
stationary state. For a broad class of single-population 
stochastic systems this approximation, complemented by 
additional perturbation techniques, yields accurate and 
controlled analytical results for the population extinction 
rate [5|. Stochastic systems involving multiple popula- 
tions present a much harder problem. Here one arrives, 
already in the leading order of the WKB approximation, 
at a generally non-integrable multi-dimensional problem 
of classical mechanics. Although the extinction rates and 
most probable paths of the system to extinction can be 
found numerically, analytical insight is usually limited. 
The situation improves, however, if the multi-population 
system exhibits time-scale separation. This may happen 
in two cases: (i) when the multi-population system is 



close to a bifurcation of the underlying deterministic rate 
equations 6, 7], and (ii) when there is a wide difference 
in individual process rates. The present paper deals with 
the second case. The process rate disparity introduces 
an additional small parameter e <C 1 which enables one 
to separate, with controlled accuracy, a two-population 
system into a fast and slow subsystems. Each of these 
subsystems is one-dimensional and therefore amenable to 
analytical solution. 

Time scale separation was recently employed in Ref. 
!S] for calculating the extinction rate of a biologically 
important component regulated by chemical reactions 
in a living cell. In that class of systems the extinction 
probability flux sets in on the slow time scale, whereas 
the fast subsystem (which rapidly adjusts to the slowly 
varying distribution of the slow subsystem) modifies the 
effective production rate of the population on the way 
to extinction. In the present work we consider a dif- 
ferent class of systems which enables us to generalize 
the standard notion of the quasi-stationary extinction 
rate by defining and calculating a time-resolved quasi- 
stationary extinction rate. It turns out that this quantity 
smoothly changes in time from a short-time asymptote 
Wi to a long-time asymptote Wi. The short-time quasi- 
stationary extinction rate W\ sets in already on the fast 
time scale. Notably, W\ coincides with the extinction 
rate in the case when the slow processes are absent al- 
together. Then, as the probability distribution evolves 
towards the long-time quasi-stationary distribution, the 
extinction rate undergoes a smooth but exponentially 
large change and approaches Wi- This evolution occurs 
on the time scale of the slow subsystem (which is much 
longer than the time scale of the fast subsystem but much 
shorter than the mean time to extinction r ex ~ W^ 1 )- 
The exponentially large disparity of W\ and W2 is an in- 
stance of the recently discovered extinction rate fragility 
, where r ex experiences a discontinuity when the rates 
of the slow processes are taken to zero. Essentially, the 
present work renders an alternative, time-resolved, de- 
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scription to the phenomenon of extinction rate fragility. 

We will present the theory on the example of a 
well-known model of epidemiology: the stochastic SIS 
(Susceptible-Infected-Susceptible) model 0, E3| with a 
slow population turnover - renewal and removal. Sec- 
tion|TT]presents the governing equations and explains how 
to exploit the time scale separation. In Section IIIII we 
approximately solve the dynamics of the fast and slow 
subsystems, calculate the time-resolved disease extinc- 
tion rate and compare our analytical predictions with a 
numerical solution of the exact master equation. Sec- 
tion [III] also deals with the mean time to extinction of a 
population exhibiting a time-dependent extinction rate. 
Section ITVl presents a discussion of our results. 



II. GOVERNING EQUATIONS AND TIME 
SCALE SEPARATION 

The stochastic SIS model 0, [l(| describes a 
Markov process involving susceptible and infected sub- 
populations. Upon recovery the infected individuals be- 
come susceptible again. The probability P n , m {t) to ob- 
serve, at time t, n susceptible and to infected individuals 
is governed by the master equation with transition rates 
from Table 1. One can always represent the renewal rate 
of susceptible individuals — an independent parameter of 
the model — as fiN, where N, an alternative independent 
parameter, scales as a typical average total population 
size in a steady state. Rescaling time by the recovery 



Event 


Type of transition 


Rate 


Infection 


S->S-l,I->I+l 


(P/N)SI 


Recovery 


S -> S + l, I -> I -1 


kI 


Renewal of susceptible 


S -> S + 1 


fj,N 


Removal of susceptible 


S->S-1 


/j.S 


Removal of infected 


J ->•/-! 


ill 



however, for pre-exponential factors, this condition can 
be relaxed to a much less restrictive one, e -C 1, as we 
explain below. 

The ultimate state of the SIS model is infection- free. 
When R — 1 = 0(1) > there is a quasi-stationary en- 
demic state with a life time T ex which is exponentially 
long with respect to N. To get an insight into how this 
quasi-stationary state is approached, let us consider the 
deterministic rate equations of the SIS model which op- 
erate with the average numbers h and to of susceptible 
and infected individuals, respectively: 



R__ 

m = — nm — m — em . 
N 



(2) 
(3) 



These equations accurately describe the dynamics of n 
and to at times short compared with r ex . The attracting 
fixed point of Eqs. ([2j and 



N(l + e) ^ N 

R ~ i? ' 



(4) 



TABLE I: Stochastic SIS model with population turnover 
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describes the endemic state of the population which is 
typically reached on the long, demographic time scale 
r e = ^> 1. (To remind the reader, our time is rescaled 
by the recovery rate k.) What happens on a much shorter 
epidemic time scale n ~ 1? It will be assumed that the 
initial average size of the total population, k = m + h, is 
greater than N / R, which is important for the further dis- 
cussion. On the time scale n ~ 1 the terms proportional 
to e can be neglected, and one can see that the average 
number of susceptible individuals approaches N/R which 
is close to . In its turn, the average number of infected 
rapidly adjusts to the current value of the average total 
population size k: rhj. = k — N/R. Restoring the e-terms, 
one can see that k = k(t) is a slow function of time, and 
its dynamics is described by the simple equation 



rate, nt — > t, one can write down the master equation for 
the SIS model as 

~T7-Rn. m(^) £-^-?n— l,m &NP n m 
at 

+e (n + 1) P n +i, m - enP ntin 

e(m + l)P, hm+1 - emP n>m 

+{R/N)(n + l){m-l)P n+l>m _ 1 - {R/N)nmP n>m 

+ (m + 1) P n -i, m+ i - mP n . m , (1) 

where e — fi/n, R — /3/k, and we assume that Pn,m — 
if at least one of the indices n or m is negative. We will 
assume throughout this work that N 3> 1 and R — 1 = 
0(1) > 0. A slow population turnover implies smallness 
of e. For our theory to be accurate it is necessary that 
a strong inequality e <C 1/N holds. If one does not care, 



k = e{N-k) 



(6) 



obtained by summing up Eqs. ^ and As a result, 
k(t) slowly flows to N, 

k{t) =N+{k ~ N)e~ et , 

and the average number of infected rhj, approaches to* 
from Eq. ([5]). On the phase plane n, fh, see Fig. [TJ the 
trajectory first rapidly reaches, the vertical line h — N/R 
and then slowly, on the long time scale r e , approaches the 
ultimate fixed point ([J} and ([5]) along this vertical line. 

How is this dctcrminisitic dynamics reflected in the 
actual evolution of P n .m(t)^ At times t\ < t < r e the 
distribution P n , m (t) is peaked at n ~ N/R,m ~ ffi^t) 
and evolves in time on the demographic time scale r e 3> 
T\. At longer times, r e < t < r ex the distribution reaches 
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FIG. 1: The phase portrait of the deterministic rate equations 
||SJ| and ([3]). The filled circles show the attracting fixed point, 
see Eqs. Q and ((5j) , and the repelling fixed point n = N,m = 
0. The arrows show the flow directions on the phase plane. 
The parameters are R — 2 and e = 5 x 10 . 



its long-time asymptote. It is important that a well- 
defined extinction probability flux sets in quite rapidly, 
at t > Ti, and it varies in time on the time scale t £ . It 
is obvious that, in general, the disease extinction rate at 
t Tg should be very different from its value at t 
r e . Indeed, even the average numbers of susceptible and 
infected are generally very different at the earlier and 
later times. It turns out, however, that an exponentially 
large disparity in the disease extinction rate at short and 
long times is observed even in the special case of k(t = 
0) = N, when the average numbers of susceptible and 
infected stay (almost) constant on the slow time scale t £ . 
The time-dependent disease extinction rate is defined 



n,m>0 



Pn 



E 



(7) 



\n,m>0 



The limit of W(t,e) as t — > oo is denoted by W 2 — T ex . 
For t < r ex one has J2 n .m>o Pn . m — As a result ) f° r 
these times 



W(t,e) 



dt 



E P n 



w(t, e), 



(8) 



n,m>0 



where w(t, e) is the disease extinction probability current. 
To calculate w(t, e) we return to the master equation 
(fTJ) and notice that the disease can disappear from the 
population only via transitions from any of the states 
(n, 1), where n = 0, 1, . . ., to the state (n, 0). Therefore, 



w(t, e) 



E 

fc=0 



U'k 



: (t) = J2(l + £)Pk-i,i(t) (9) 



k=l 



where 



w k (t) 




ePjfe.i, k > 0, 
k = 0. 



(10) 



As a result, for t <C r ex 

oo 

W{t, e) ~ w(t, e) = 5^(1 + e) P fc _ x , i(t). (11) 
fe=i 

Before exploiting the time-scale separation, let us 
first obtain some exact relations. Consider Pk = 
^ j Pk-m, m (k > 1) which is the probability to find 
k > 1 individuals so that at least one of them is infected. 
Summing Eq. (fTJ) over m while keeping the total number 
k = n + to of individuals constant, we obtain an exact 
equation 



dt 



P k {t) = eNP k _ 1 -e(N + k)P k +£(k + l)P k+1 



W k -l 



(12) 



where Pq — 0. An additional exact equation can be ob- 
tained once we represent the probability P n ,m{t) as 



m ~~ n m _[_ n _ £ m >i 



|m 5 



(13) 



where n m+n= /- im >i \ m (t) = P n , m {t) [Pfe(*)] , defined 
for m > 1, is the probability to have to infected individ- 
uals conditioned on n + to = k and m > 1. We will use 
the following shorthand: 



Pkijn) — ^-m+n=k,m>l \m(t) 



(14) 



This conditional probability is identically equal to for 
m < and to > k, and it obeys the following exact 
equation: 

^Pfe(m) = P k {\) P k {m) + eP k+l {\) P k {m)P k+1 (P^ 1 
R 

+ — (k -m+ l)(m - l)Pfc(m - 1) 



(fc — m)mP k (m) — mP k (m) + (m + 1) P k (m + 1) 



i? 

- (P fc ) _1 j^Pfc-i [pcM - Pfe-i(m) 
+e (fc + 1) P k+1 \p k {m) - P fc+ i(m 

+eP k +i 



mP k+1 (m) - (to + l)P fc+ i(TO + 1) 
The disease extinction rate (|11[) in the new notation reads 

oo 

W{t, e) ~ ttf(t, e) = (1 + e) ^ P fc (t)P k (l)(t) . (16) 



(15) 



fc=i 



At t r ex , the exponentially small term u>fc-i in 
Eqs. (fT2"f can be neglected, and we obtain 



fl fc (t) = eJVflb_i-e(JV + A)flfc + e(fe + l)fl t+ i. (17) 



We see that the evolution of P k (t) at i <§; r ex proceeds on 
the slow time-scale r e and is decoupled from the evolution 
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of Pfe(m). Equation (| 1 T|) describes the relaxation of the 
probability distribution P k (t) to the steady state distri- 
bution . The latter is determined from the equation 



eNP, 



(o) 
fc-i 



£ (N + k)P^+e(k + l)Pj;%=0 : 



(18) 



subject to normalization X^feLi Rk°^ = 1 
The evolution of Pk(m) is fast at t « 
ti < t < r ea; . The solution of Eq. {Tf[ 
can be sought in the form 



7~i and slow at 
at n « i < Tex 



P k (m,t) = pj?\m) + eNPj; 1 \m,et) 
+ {eN) 2 Pt 2) (m,st) + ... . 

where P^ (m) obeys the stationary equation 
^ (k - m + l)(m - 1) Pf } (m - 1) 



(19) 



Pfe 

— {K — 111)111 1 k \„ t j — ,,v± k 

+ (m + l)P^>(m + 1) + P fe w (l)P fc w (m) = 0, (20) 



[k — m)m Pj,°\m) — m Pj,°\m) 



P 

=i ^fc 



(0) 



where 

can confine ourselves only to the leading term P^°\m) 



(m) = 1. For eN <C 1 and f > n we 



in Eq. (I19[) . In the following we will omit the superscript 
(0). Importantly, Eq. (|2"0"j) does not include P k . 

Each of the decoupled equations (fTTf and (|20p has 
a simple meaning. Equation (|20[) describes a one- 
dimensional quasi- stationary distribution (QSD) of the 
number of infected in the SIS model without population 
turnover, where the total population size is conserved, 
k = const. This QSD forms relatively rapidly, at t > n, 
so we can call this subsystem fast. Once the QSD is 
found, Pfe(l) yields the disease extinction rate for the 
given k. 

In its turn, Eq. (fP7|) describes the evolution of a one- 
dimensional time- dependent distribution of the total pop- 
ulation size k. The characteristic time scale of this time 
dependence is t £ = 3> Ti, so we can call this sub- 
system slow. Having found the slowly evolving distribu- 
tion, one can calculate the time-resolved disease extinc- 
tion rate from Eq. (fT6|) (with the e-term dropped to avoid 
excess of accuracy): 



w(t,s)~^Tp k (t)p k (i) 



(21) 



k=l 



This rate is nothing but the average of the instantaneous 
extinction rate for a given k (found from the fast subsys- 
tem) over the time-dependent fc-distribution (found from 
the slow subsystem). This result is valid when f > ti 
and eN < 1. 

The next section deals with the solution of Eqs. (fT7|) 
and (|20|) . and with calculating W(t, e), for a class of ini- 
tial conditions for which the total number of individuals 
at t — is equal to N: 



5k. 



N 



(22) 



where 6 k ,N is the Kronecker's delta. This special choice of 
initial condition will make it possible to relate the fore- 
going time-dependent picture of extinction to the phe- 
nomenon of extinction rate fragility. 



III. SOLVING THE TIME-SCALE-SEPARATED 
PROBLEM 

A. Fast subsystem 

Disease extinction in the stochastic SIS model without 
population turnover has been extensively studied starting 
from the pioneering paper by Weiss and Dishon 10]. The 
mean time to extinction of the disease in this case was 
first obtained by Nasell [ll|, see also Refs. [1, [l2j- The 
extinction rate, rescaled to the recovery rate K, is equal 
to 



P fc (l) 



* (Rk - 1) 5 



2vr R k 



x exp 



1 

Rk 



lnP fc -l 



(23) 



where R k = kR/N; here and in the following the su- 
perscript (0) in (1) is omitted. Equation (|2"3"]l holds 
when the factor in the exponent is sufficiently large in 
absolute value 0, 0j], (3 : 



1 

Rk 



lnP fc 



1 > 1. 



The particular value of the extinction rate for k — N is 
nothing else but W\ : the disease extinction rate observed 
at times T\ < t <C t £ , when the distribution of infected 
has already adapted to the current value of k, but the 
fc-distribution has not yet evolved and is still close to the 
Kronecker's delta (l22l: 



Wi 



N (R - 1) 



2tt R 



x exp 



V [ -+hiP-l 



(24) 



B. Slow subsystem 



Equation (TiT]) is exactly solvable 13( with the help of 
the probability generating function 



G(z,T)=J2z k Pk(T). 



(25) 



k=0 



where z is an auxiliary variable, and t = t/r e = st. Once 
G(z, t) is found, the probabilities Pfe(r) can be recovered 
from the Taylor expansion: 



Pk(r) 



1 d k G(z,r) 
k\ dz k 



(26) 
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After a simple algebra the master equation (ITTl) becomes 
an evolution equation for G(z,r): 



dG 
dr 



ti-.)(i-™ 



(27) 



This first-order PDE can be solved by characteristics. 
The general solution is 



G(z,r) = e Nz f 



1 - z 



(28) 



where /(£) is an arbitrary function. To determine /(£), 
one should use the initial condition. In terms of Pk it is 
given by Eq. (|2"2")l . In terms of G we obtain 



G(z,t = 0) = 



iV 



This yields /(£) = (1 -^^e-^ 1 ^, and so the resulting 
solution, in terms of G, is 



G(z,r) 



1 



1 



N 



exp[7V(z- 1)(1 



(29) 

In the limit of t > 1 we obtain G(z, t) = exp[N(z — 1)]. 
This describes a Poisson distribution with mean N: 



Pfe(Or £ 



1 d k G(z lT = oo) 



dz k 



N k e 



-N 



z=0 



k\ 



(30) 



Using this Poisson distribution we can calculate Wz 



-l. 



the long-time asymptote of the disease extinction 
rate W(t e <C i < T el ), observed when the fc-distribution 
has already reached quasi-stationarity. As k 3> 1, we can 
use Stirling's formula for fc! and obtain 



P fe (r e < t < r ex ) 



^N(x—l—x\nx) 



V2ttNx 



(31) 



where x = k/N . Now we calculate W(t ^> r e ) = W 2 by 
using the distributions (f2"3"|) and (j3"Tj) . Replacing the sum 
over k by an integral over x — k/N in Eq. (|16p. we obtain 



Wo 



N f , {Rx-lf 

— I ax 

2tt J Rx 

R- 1 



x exp [-N {x In Rx + x In x - 2x + 1 + 1 / R)] . (32) 

As A ^> 1, the integral can be evaluated by the saddle- 
point method (as a result, the exact location of the lower 
bound of integration is actually unimportant). The sad- 
dle point condition \n(Rx s ) +ln x s = yields x s = 1/ V^R, 
corresponding to fc = N/y/R. By virtue of Eq. ([2"T]) . 
/j = N/yf~R is the most probable total population size 
when the number of infected is exactly one. Performing 
the gaussian integration, we obtain 



W 2 



N(yjR-l) 2 
2 v ^Fi? 3 / 4 



exp 



-N 1 



1 



(33) 



This result, without the pre-exponential factor, was ob- 
tained by Khasin and Dykman [9j who used (the lead- 
ing order of) WKB approximation. The important pre- 
exponent has not been known previously. 

For the special class of initial conditions (|22|) . the 
short-time, Wi, Eq.(g3]), and the long-time, W 2 , Eq. (j3"3"l) . 
quasistationary extinction rates correspond to the qua- 
sistationary extinction rates observed without population 
turnover and with a very slow population turnover, re- 
spectively. For R — 1 = O(l) the rate W 2 is exponen- 
tially larger than W\ . This exponential disparity reflects 
fragility of the extinction rate of the system without pop- 
ulation turnover with respect to addition of slow popula- 
tion turnover. Although the fragility phenomenon was 
established in the e-domain Q, we see that a closely 
related phenomenon is also observed, for proper initial 
conditions, in the time domain. 

Now we discuss an important applicability criterion of 
our theory. The main contribution to the integral (|32[) 
comes from a relatively narrow, 0(y/N), vicinity of the 
saddle point k = N/ \J~R~. In this vicinity, Eq. (|23|) is valid 
if it describes a quasi-stationary distribution on the time 
scale t s . The corresponding criterion, e 3> Pjq/^/^X) 
serves as the lower bound on e for the validity of Eq. (|33l) . 
For N > 1 and R - 1 = 0(1) the quantity P N/ ^{1) is 
exponentially larger than W 2 . Therefore, the criterion 
£ S> Pjsjj^/r^X) i s much more restrictive than the obvious 
criterion e 3> W 2 . Returning to the fragility problem, 
we note that at e ^ -P/v/,/r(1) the extinction rate experi- 
ences a gradual crossover from W 2 to W\ in the e-domain. 

Let us return to Eq. (j2"9")) and find the time-dependent 
distribution i\ (r) from Eq. (|26| . It is convenient to cal- 
culate the derivatives in the complex z-plane by using 
the Cauchy theorem: 



— G(z,t), 



(34) 



where the integration contour encircles the pole z — 
of the complex plane z. As k ^> 1 and iV 3> 1, we can 
evaluate the contour integral using the saddle-point ap- 
proximation and deforming the contour so as to achieve 
the steepest descent. Using Eq. (|2U)) , we rewrite Eq. (|3"4")l 
as 



Pk(r) 



1 

2vri 



dz ■ 



N&(z,x,t) 



(35) 



where 



z-l 



$ = In \ l + — —j +(z- 1)(1 - e~ T ) -xlnz 

and x = k/N. The saddle point z = z*(x,t) is de- 
termined from the equation d^/dz = 0. The suitable 
solution is 

z»(x,r) = [2(l-e- r )]- 1 [ 1 + 2 ; -2coshr 



y/3 + x(x - 6) + 4(x - 1) cosh r + 2 cosh 2t .(36) 
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Now we expand 



1 



At t > one has z, r) > 0. Therefore, the steep- 

est descent of function $ occurs along the straight line 
Kez — of the complex plane, see Fig. [5J so we de- 
form the contour accordingly. In view of ./V ^> 1 a small 
segment of the straight line Re z — z„ gives a dominant 
contribution to the integral, and the gaussian integration 
yields 

p N$(z, ,x,t) 

Pk(r) , . (37) 

As a simple test of this result, one can go to the limit 
r ^> 1 and obtain and $" = 1/x. Then Eq. (|37| 

yields Eq. (T5Tj) as expected. 




Rez 



FIG. 2: The steepest descent path for Eq. 



Time-resolved extinction rate 



Now we return to Eq. (12ip and calculate the time- 
resolved extinction rate W(t, e) = W(t) of the disease 
by averaging the instantaneous extinction rate (|23[) over 
the slowly time-dependent distribution (|3"T)) . Replacing 
the sum over k by an integral over x = k/N and evalu- 
ating the integral by the saddle-point method, we obtain 
after some algebra: 



W{t) 



[Rx s (i 
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Nx s (t) 



Rx s (t)z4x s (t),t] y 2it<I>"[x s (t),t}A"[x s (t),t} 

x e -N[R- 1 +x s (T)(lnx s (r)B.~l)+S s ^Ar),r)} _ /ggs 

Here the time-dependent saddle point .t = x s (t) is de- 
termined by the equation 



x s R [l + c(x s ,T)e T ] = 1 



(39) 



where 



2 sinh r 



c(x,r)=-[2(e T -l)]- 1 

-Vtl - x + 2 sinh r) 2 - 4(1 - x) (e T - 1) . (40) 
Furthermore, the function S s (x,t) is given by 
S s (x,t) = (1 + c)" 1 [(1 + c-c 2 ) In (l + ce T ) 

+ c(l + c)e T [In (1 + ce T ) - 1] - ce~ T In (1 + ce T ) 
+ c-ln(l + c) (41) 



with c = c(x,t) from Eq. ((40)) . Finally, 

A(jg,t) = iT 1 +x(lnxi?- 1) + S s {x, r) , 

and we have written for brevity [z* (x s , r) , x s (r) , r] = 
$"Mt),t]. 

Equation (|38| is one of the main results of this work. 
It describes a time- resolved disease extinction rate Wist) 
which smoothly changes from W\ at t <£L r e to TT2 at t ^ 
r e . Figures [3] and d] show, by solid lines, typical examples 
of this behavior for two different sets of parameters. 

To test the theoretical time-dependent extinction rate 
W(et) predicted by Eq. (f3"5)) we solved numerically a 
(truncated version of the) original master equation ((T|), 
using the ODE45 solver of MATLAB. A very good agree- 
ment between the theory and the numerical solution 
of the original master equation ([T]) was obtained for 
TV = 200 and e = 10~ 3 , see Fig. [3l Here the time scale 
separation criterion eN -C 1 was satisfied. 

Importantly, the (quite restrictive) criterion eN <C 1 
can be replaced by a much less restrictive one e <C 1 if one 
does not care about pre-exponential factors in Eq. (J38J). 
Indeed, criterion eN -C 1 appears when one implements 
the time-scale-separation procedure directly in the mas- 
ter equation |T]). One can follow a different strategy, 
however, and start with applying a time-dependent WKB 
approximation to the full two-dimensional master equa- 
tion ([TJ) . A proper WKB ansatz is 

= a{x,y,t) exp[-NS(x,y,t)], (42) 

where x = n/N and y — m/N. In the leading WKB- 
order one obtains a two-dimensional time-dependent 
Hamilton- Jacobi equation. The corresponding Hamilto- 
nian H(x,y,p x ,p y ,t;e) is independent of N. The small 
parameter e<1, present in the Hamiltonian, allows one 
to perform a time-scale-separation procedure by seeking, 
for t > Ti, 

S(x,y,t) = S (x 1 y,et)+eSi{x 1 y,et)+e 2 S 2 {x,y,et) + . . . , 

(43) 

where So ~ Si ~ . . . ~ 1, and So(x,y,et) has a separa- 
ble structure. Our derivation, leading to Eq. (|3"5)) . yields 
So(x,y,et) and a{x,y,et) but not Si, 52, The eSi 
term makes a contribution of order eN in the exponent 
of Eq. (|32]l. This contribution is negligible if eN « 1. 
In this case the pre-exponent in Eq. (f3"5| is accurate as 
we have already seen. On the contrary, if eN > 1, the 
unknown correction eS\{x, y, et) becomes significant, and 
the account of the pre-exponent in Eq. (|38p would be in 
excess of accuracy. Now, what happens if we only need 
an accurate estimate of lnW(t,e)/N at JV 3> 1? Here 
both the pre-exponent in Eq. (|38[) . and the unknown cor- 
rection eSi (x, y, et) become negligible, and the result is 
described by the exponent in Eq. (|3"51) . Indeed, we ob- 
served an excellent agreement between theory and nu- 
merical calculations for \nW(t,e) when the parameter 
eN was comparable with, or even larger than 1, see for 
example Fig. 2J The numerical results, presented in Figs. 
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[3] and [4] are converged both with respect to the trunca- 
tion of the master equation, and with respect to the error 
tolerance of the MATLAB solver. 




FIG. 3: The time-resolved disease extinction rate W versus 
the rescaled time ent = fit for N — 200, R — 10 and e — 



10 as predicted by Eq. (|38[) (solid line) and obtained by 
a numerical solution of the (truncated) master equation JT|) 
(the dashed line). 




FIG. 4: The natural logarithm of the the time-resolved disease 
extinction rate W versus the rescaled time ent = fit for N — 
100, R = 4 and e = 10~ 2 as predicted by Eq. (38} (solid 
line) and obtained by a numerical solution of the (truncated) 
master equation (TjTJ (the dashed line). 



D. Mean time to extinction 

The time-resolved extinction rate, which we have cal- 
culated in this work, provides a sharp characterization 
of the stochastic population dynamics. This character- 
ization is lost if one is only interested in the average 
extinction quantities such as the mean time to extinc- 
tion (MTE) r ex of the population. To better understand 
this point, consider the disease extinction probability as 
a function of time: 



At times, t% < t <C r ex , the growth rate of of Vo(t) obeys 
the relation 



dV (t) 
dt 



W(t,e) 



(44) 



which follows from Eq. (|S} and the conservation of the to- 
tal probability. At these times the extinction probability 
rate experiences a smooth but exponentially large change 
with time on the time scale of r E . 

At longer times, t ~ r ex » r e , Eq. (|44l) no longer holds, 
and should be replaced by the relation 



dV (t) 
dt 



W 2 e 



-W 2 t 



(45) 



The mean time to extinction r ex can be obtained by av- 
eraging t over dPo(t)/dt (which is the probability distri- 
bution of extinction times) : 



t^dt. 
dt 



(46) 



The dominant contribution to this integral comes from 
times t 3> r e where <fP (i) / dt is determined by Eq. (|45p . 
Therefore, up to exponentially small corrections, 



r W 2 te- W2t dt 
Jo 



1/W 2 . 



(47) 



That is, the time-resolved extinction rate Wit, e) pro- 
vides a much more detailed information about the system 
dynamics than the MTE (which only probes the late-time 
asymptote of the extinction rate). 



IV. DISCUSSION 

We have addressed extinction of a population in a 
two-population system in the case when the population 
turnover - renewal and removal - is much slower than all 
other processes. The ensuing time scale separation makes 
it possible to introduce a short-time quasi-stationary ex- 
tinction rate W\ and a long-time quasi-stationary extinc- 
tion rate W 2 , and develop a time-dependent theory of 
the smooth transition between the two rates. The quan- 
tities W\ and W 2 coincide with the extinction rates when 
the population turnover is absent altogether, and present 
but very slow, respectively. The exponential difference 
between the two rates reflects fragility of the extinction 
rate in the population dynamics without turnover 
The present work suggests an alternative picture of the 
extinction rate fragility: in the time domain instead of 
the e-domain where it was originally established. 

Our main results can be expressed in the following way. 
We started out by presenting the probability distribution 
P n ,m>o(t) in a factorized form: P n , m >o(t) = Pk(t)P k {m), 
where Pfc(t) is the probability to have the total popula- 
tion size k, when at least one individual is infected, and 
Pfe(m) is the probability to have m > infected indi- 
viduals, when the total population size is k. At t 3> 71 
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and eN <C 1 the Pk (m)-distribution is independent of 
time to the zero order in sN; the population turnover 
only affects Pk(t). As a result, the time-dependent ex- 
tinction rate W(t, e) is determined by the extinction rate 
for the total population size fc, obtained for e = and av- 
eraged over the slow-time-dependent Pfc (i)-distribution. 
The short-time quasi-stationary extinction rate W\ corre- 
sponds to the initial probability distribution of the total 
population size fc, whereas the long-time quasi-stationary 
extinction rate W 2 corresponds to the steady-state k- 
distribution. Under less restrictive conditions e <C 1 and 
N ^S> 1 our theory accurately predicts the logarithm of 
the time-dependent extinction rate hiW(t,e). 

We have shown that the time-resolved quasi-stationary 
extinction rate encodes a more detailed information 
about the stochastic dynamics than the average quan- 
tities such as the mean time to extinction. The latter 
quantity is determined by the late-time asymptote W2 of 
the time-resolved extinction rate. 

Our analytical approach can be used in a host of 
two-population models which exhibit a long-lived quasi- 
stationary state on the way to extinction, and where a 
disparity between the process rates enables one to sep- 
arate the system into two one-dimensional sub-systems: 
the fast and the slow. One can envision two different sce- 
narios in such systems. In one scenario, extinction takes 
place only in the slow subsystem, whereas the fast sub- 
system merely modifies the effective process rates, as in 



Ref. Q. In another scenario the problem is reducible (as 
in the SIS model with a slow population turnover which 
we have considered here) to averaging the instantaneous 
extinction rate, generated by the fast subsystem, over 
the time-dependent distribution of the slow subsystem. 
For a sufficiently large population size, the instantaneous 
extinction rate, generated by the fast subsystem, can be 
accurately calculated via WKB approximation [5j. How 
can one find the time-dependent population size distribu- 
tion of the slow sub-system? For the SIS model the one- 
dimensional master equation (1171) for the slow population 
turnover is decoupled from the fast sub-system, and even 
exactly solvable. In general one cannot hope for such a 
dramatic simplification, and the fast subsystem will mod- 
ify the effective process rates of the slow subsystem, as 
in Ref. 8]. It is important, however, that this effect can 
be described by a time-dependent WKB theory, once the 
corresponding large parameter is present. 
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